install.packages("MASS")
library(MASS)
install.packages("ggplot2")
library(ggplot2)


#simulation based on data on period #3 (1999-2021)
sacL.resc1999 <- sacsarlm(lnme_abs ~ lnme_abslag + lngdp_abs + lnpop + as.factor(ccode)+counter+
                            counter.sq+counter.cub, data = data.ordered1999, capdistinv.2listw.nors.resc1999, Durbin = F) #weights non-normalized
sacL.resc1999$rho

# VCOV (matrix)
vcov <- vcov(sacL.resc1999)
#View(vcov)

VCVM <- vcov[-c(2:3,5:35),-c(2:3,5:35)] #let's select relevant parameters: rho and phi 
VCVM # rho, phi


draws <- 1000
coef(sacL.resc1999)[1] # rho
coef(sacL.resc1999)[4] # phi
simdraws <- mvrnorm(draws,c(coef(sacL.resc1999)[1],coef(sacL.resc1999)[4]),VCVM) #library(MASS)
#View(simdraws)

fx_boot <- matrix(NA, nrow=draws, ncol=3) #to contain simulated average long-run effects (direct, total, and indirect)
direct_avg_LR_boot <- matrix(NA, nrow=18, 1) #to contain cumulative direct effects, over the period, for 18 separate countries  

#[1] "United States"  "Canada"         "United Kingdom" "Netherlands"    "Belgium"        "Luxembourg"     "France"        
#[8] "Spain"          "Portugal"       "Germany"        "Poland"         "Hungary"        "Czech Republic" "Italy"         
#[15] "Greece"         "Norway"         "Denmark"        "Turkey" 

sims.us = matrix(NA, nrow=draws,ncol=23)
sims.others1=matrix(NA,nrow=draws, ncol=(12*23)) #multiplier 12 in (12*23) because from 2 to 13, there are 12 countries (each having 23 year periods)
sims.it = matrix(NA, nrow=draws,ncol=23)
sims.el= matrix(NA, nrow=draws,ncol=23)
sims.others2=matrix(NA, draws, (2*23)) #countries #16 and #17
sims.tur = matrix(NA, nrow=draws,ncol=23)

for(k in 1:draws){
  estcf1 <- (diag(566) - simdraws[k,2]*TL - simdraws[k,1]*capdistinv.wide1999.resc)  #I - phi*L - rho*W (based on Cook et al. 2023)
                                                                                     ##TL is lag multiplier (L), 1999-2021
  estcf <- solve(estcf1) #inv(I - phi*L - rho*W); we simulate 1000 spatial effects matrices 
  
  for(j in 1:18){ #we here collect in matrix 'direct_avg_LR_boot' cumulative direct effects, over the period, for 18 separate countries  
    #direct_avg_LR_boot[j,] <- sum(estcf[,j][seq(j, nt, n)] ) 
    if(j < 14){
      #[1] "United States"  "Canada"         "United Kingdom" "Netherlands"    "Belgium"        "Luxembourg"     "France"        
      #[8] "Spain"          "Portugal"       "Germany"        "Poland"         "Hungary"        "Czech Republic"  
      
      direct_avg_LR_boot[j,] <- sum(estcf[,j][c(seq(j, 90, 18),seq((90+j),215,25),seq((215+j), 566,27))]) 
    } else {
      if(j == 14){ #"Italy"  
        direct_avg_LR_boot[j,] <- sum(estcf[,j][c(seq(j,90,18),seq((91+j),215,25),seq((216+j), 566,27))])
      }else{
        if(j == 15){ #"Greece" 
          direct_avg_LR_boot[j,] <- sum(estcf[,j][c(seq(j,90,18),seq((92+j),215,25),seq((219+j), 566,27))])
        }else{ # "Norway"         "Denmark"        "Turkey" 
          direct_avg_LR_boot[j,] <- sum(estcf[,j][c(seq(j,90,18),seq((97+j),215,25),seq((224+j), 566,27))])
        }}}                
  }#loop within a loop
  
  fx_boot[k,1] <- mean(direct_avg_LR_boot) #average direct long run effect
  fx_boot[k,2] <- mean(colSums(estcf[,1:18])) #average long run total effect
  fx_boot[k,3] <- fx_boot[k,2] - fx_boot[k,1] #average long run indirect
  
  for(i in 1:18){ #in country-specific matrices, we collect indirect effects (a given country has on other allies) 
                  ##for each year throughout the entire period from 1999 to 2021, hence 23 year-periods
    if(i==1){
      sims.us[k,]=  cbind(sum(estcf[2:18,i]),  sum(estcf[20:36,i]), sum(estcf[38:54,i]),  sum(estcf[56:72,i]),  sum(estcf[74:90,i]),
                          sum(estcf[92:115,i]),  sum(estcf[117:140,i]), sum(estcf[142:165,i]),   sum(estcf[167:190,i]),   sum(estcf[192:215,i]),
                          sum(estcf[217:242,i]),  sum(estcf[244:269,i]), sum(estcf[271:296,i]), sum(estcf[298:323,i]), sum(estcf[325:350,i]),
                          sum(estcf[352:377,i]), sum(estcf[379:404,i]), sum(estcf[406:431,i]), sum(estcf[433:458,i]), sum(estcf[460:485,i]),
                          sum(estcf[487:512,i]),sum(estcf[514:539,i]), sum(estcf[541:566,i]))
    }else{
      if(i<14){
        sims.others1[k,((23*(i-1)-22):(23*(i-1)))] = cbind(sum(estcf[1:(i-1),i],estcf[(i+1):18,i]), #1st year
                                                           sum(estcf[19:(19+i-2),i],estcf[(19+i):36,i]), #2nd year
                                                           sum(estcf[37:(37+i-2),i],estcf[(37+i):54,i]),
                                                           sum(estcf[55:(55+i-2),i],estcf[(55+i):72,i]),
                                                           sum(estcf[73:(73+i-2),i],estcf[(73+i):90,i]), 
                                                           sum(estcf[91:(91+i-2),i],estcf[(91+i):115,i]), 
                                                           sum(estcf[116:(116+i-2),i],estcf[(116+i):140,i]),
                                                           sum(estcf[141:(141+i-2),i],estcf[(141+i):165,i]),
                                                           sum(estcf[166:(166+i-2),i],estcf[(166+i):190,i]),
                                                           sum(estcf[191:(191+i-2),i],estcf[(191+i):215,i]),
                                                           sum(estcf[216:(216+i-2),i],estcf[(216+i):242,i]),
                                                           sum(estcf[243:(243+i-2),i],estcf[(243+i):269,i]),
                                                           sum(estcf[270:(270+i-2),i],estcf[(270+i):296,i]),
                                                           sum(estcf[297:(297+i-2),i],estcf[(297+i):323,i]),
                                                           sum(estcf[324:(324+i-2),i],estcf[(324+i):350,i]),
                                                           sum(estcf[351:(351+i-2),i],estcf[(351+i):377,i]),
                                                           sum(estcf[378:(378+i-2),i],estcf[(378+i):404,i]),
                                                           sum(estcf[405:(405+i-2),i],estcf[(405+i):431,i]),
                                                           sum(estcf[432:(432+i-2),i],estcf[(432+i):458,i]),
                                                           sum(estcf[459:(459+i-2),i],estcf[(459+i):485,i]),
                                                           sum(estcf[486:(486+i-2),i],estcf[(486+i):512,i]),
                                                           sum(estcf[513:(513+i-2),i],estcf[(513+i):539,i]),
                                                           sum(estcf[540:(540+i-2),i],estcf[(540+i):566,i]))#23rd year
      }else{
        if(i==14) {
          sims.it[k,] = cbind(sum(estcf[1:(i-1),i],estcf[(i+1):18,i]),
                              sum(estcf[19:(19+i-2),i],estcf[(19+i):36,i]),
                              sum(estcf[37:(37+i-2),i],estcf[(37+i):54,i]),
                              sum(estcf[55:(55+i-2),i],estcf[(55+i):72,i]),
                              sum(estcf[73:(73+i-2),i],estcf[(73+i):90,i]), 
                              sum(estcf[91:(91+i-1),i],estcf[(91+i+1):115,i]), 
                              sum(estcf[116:(116+i-1),i],estcf[(116+i+1):140,i]),
                              sum(estcf[141:(141+i-1),i],estcf[(141+i+1):165,i]),
                              sum(estcf[166:(166+i-1),i],estcf[(166+i+1):190,i]),
                              sum(estcf[191:(191+i-1),i],estcf[(191+i+1):215,i]),
                              sum(estcf[216:(216+i-1),i],estcf[(216+i+1):242,i]),
                              sum(estcf[243:(243+i-1),i],estcf[(243+i+1):269,i]),
                              sum(estcf[270:(270+i-1),i],estcf[(270+i+1):296,i]),
                              sum(estcf[297:(297+i-1),i],estcf[(297+i+1):323,i]),
                              sum(estcf[324:(324+i-1),i],estcf[(324+i+1):350,i]),
                              sum(estcf[351:(351+i-1),i],estcf[(351+i+1):377,i]),
                              sum(estcf[378:(378+i-1),i],estcf[(378+i+1):404,i]),
                              sum(estcf[405:(405+i-1),i],estcf[(405+i+1):431,i]),
                              sum(estcf[432:(432+i-1),i],estcf[(432+i+1):458,i]),
                              sum(estcf[459:(459+i-1),i],estcf[(459+i+1):485,i]),
                              sum(estcf[486:(486+i-1),i],estcf[(486+i+1):512,i]),
                              sum(estcf[513:(513+i-1),i],estcf[(513+i+1):539,i]),
                              sum(estcf[540:(540+i-1),i],estcf[(540+i+1):566,i]))
          
        }else{
          if(i==15){ #which(data.ordered1999$country == "Greece") 
            sims.el[k,] = cbind(sum(estcf[1:(i-1),i],estcf[(i+1):18,i]), 
                                sum(estcf[19:(19+i-2),i],estcf[(19+i):36,i]),
                                sum(estcf[37:(37+i-2),i],estcf[(37+i):54,i]),
                                sum(estcf[55:(55+i-2),i],estcf[(55+i):72,i]),
                                sum(estcf[73:(73+i-2),i],estcf[(73+i):90,i]), 
                                sum(estcf[91:(91+i),i],estcf[(91+i+2):115,i]), 
                                sum(estcf[116:(116+i),i],estcf[(116+i+2):140,i]),
                                sum(estcf[141:(141+i),i],estcf[(141+i+2):165,i]),
                                sum(estcf[166:(166+i),i],estcf[(166+i+2):190,i]),
                                sum(estcf[191:(191+i),i],estcf[(191+i+2):215,i]),
                                sum(estcf[216:(216+i+2),i],estcf[(216+i+4):242,i]),
                                sum(estcf[243:(243+i+2),i],estcf[(243+i+4):269,i]),
                                sum(estcf[270:(270+i+2),i],estcf[(270+i+4):296,i]),
                                sum(estcf[297:(297+i+2),i],estcf[(297+i+4):323,i]),
                                sum(estcf[324:(324+i+2),i],estcf[(324+i+4):350,i]),
                                sum(estcf[351:(351+i+2),i],estcf[(351+i+4):377,i]),
                                sum(estcf[378:(378+i+2),i],estcf[(378+i+4):404,i]),
                                sum(estcf[405:(405+i+2),i],estcf[(405+i+4):431,i]),
                                sum(estcf[432:(432+i+2),i],estcf[(432+i+4):458,i]),
                                sum(estcf[459:(459+i+2),i],estcf[(459+i+4):485,i]),
                                sum(estcf[486:(486+i+2),i],estcf[(486+i+4):512,i]),
                                sum(estcf[513:(513+i+2),i],estcf[(513+i+4):539,i]),
                                sum(estcf[540:(540+i+2),i],estcf[(540+i+4):566,i]))       
          }else{ 
            if(i==18){
              #which(data.ordered1999$country == "Turkey") 
              sims.tur[k,] = cbind(sum(estcf[1:(i-1),i]), 
                                   sum(estcf[19:(19+i-2),i]),
                                   sum(estcf[37:(37+i-2),i]),
                                   sum(estcf[55:(55+i-2),i]),
                                   sum(estcf[73:(73+i-2),i]), 
                                   sum(estcf[91:(91+i+5),i]), 
                                   sum(estcf[116:(116+i+5),i]),
                                   sum(estcf[141:(141+i+5),i]),
                                   sum(estcf[166:(166+i+5),i]),
                                   sum(estcf[191:(191+i+5),i]),
                                   sum(estcf[216:(216+i+7),i]), 
                                   sum(estcf[243:(243+i+7),i]),
                                   sum(estcf[270:(270+i+7),i]),
                                   sum(estcf[297:(297+i+7),i]),
                                   sum(estcf[324:(324+i+7),i]),
                                   sum(estcf[351:(351+i+7),i]),
                                   sum(estcf[378:(378+i+7),i]),
                                   sum(estcf[405:(405+i+7),i]),
                                   sum(estcf[432:(432+i+7),i]),
                                   sum(estcf[459:(459+i+7),i]),
                                   sum(estcf[486:(486+i+7),i]),
                                   sum(estcf[513:(513+i+7),i]),
                                   sum(estcf[540:(540+i+7),i]))        
            }else{
              { #e.g., which(data.ordered1999$country == "Denmark") 
                sims.others2[k,((23*(i-15)-22):(23*(i-15)))] = cbind(sum(estcf[1:(i-1),i],estcf[(i+1):18,i]), 
                                                                     sum(estcf[19:(19+i-2),i],estcf[(19+i):36,i]),
                                                                     sum(estcf[37:(37+i-2),i],estcf[(37+i):54,i]),
                                                                     sum(estcf[55:(55+i-2),i],estcf[(55+i):72,i]),
                                                                     sum(estcf[73:(73+i-2),i],estcf[(73+i):90,i]), 
                                                                     sum(estcf[91:(91+i+5),i],estcf[(91+i+7):115,i]), 
                                                                     sum(estcf[116:(116+i+5),i],estcf[(116+i+7):140,i]),
                                                                     sum(estcf[141:(141+i+5),i],estcf[(141+i+7):165,i]),
                                                                     sum(estcf[166:(166+i+5),i],estcf[(166+i+7):190,i]),
                                                                     sum(estcf[191:(191+i+5),i],estcf[(191+i+7):215,i]),
                                                                     sum(estcf[216:(216+i+7),i],estcf[(216+i+9):242,i]),
                                                                     sum(estcf[243:(243+i+7),i],estcf[(243+i+9):269,i]),
                                                                     sum(estcf[270:(270+i+7),i],estcf[(270+i+9):296,i]),
                                                                     sum(estcf[297:(297+i+7),i],estcf[(297+i+9):323,i]),
                                                                     sum(estcf[324:(324+i+7),i],estcf[(324+i+9):350,i]),
                                                                     sum(estcf[351:(351+i+7),i],estcf[(351+i+9):377,i]),
                                                                     sum(estcf[378:(378+i+7),i],estcf[(378+i+9):404,i]),
                                                                     sum(estcf[405:(405+i+7),i],estcf[(405+i+9):431,i]),
                                                                     sum(estcf[432:(432+i+7),i],estcf[(432+i+9):458,i]),
                                                                     sum(estcf[459:(459+i+7),i],estcf[(459+i+9):485,i]),
                                                                     sum(estcf[486:(486+i+7),i],estcf[(486+i+9):512,i]),
                                                                     sum(estcf[513:(513+i+7),i],estcf[(513+i+9):539,i]),
                                                                     sum(estcf[540:(540+i+7),i],estcf[(540+i+9):566,i]))   
              } }}}}}
    
  }
}

#let's collect the estimates into one place
all.countries.means = matrix(NA, 23, 18) #23 periods and 18 allies (as of 1999)
countries.list = as.vector(c("US","CAN", "UK","NL","BE","LU","FR","ES","PT","DE","PL","HU","CZ","IT","EL","NO","DK","TUR"))
colnames(all.countries.means) = countries.list

all.countries.upr = matrix(NA, 23, 18) #matrix to contain upper levels of 95% CI
colnames(all.countries.upr) = countries.list

all.countries.lwr = matrix(NA, 23, 18) ##matrix to contain lower levels of 95% CI
colnames(all.countries.lwr) = countries.list

#
#US
#
us.cum =matrix(NA,23, draws) #to contain cumulative (over time) indirect effects
t.sims.us=t(sims.us)
dim(t.sims.us)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      us.cum[j,i] = t.sims.us[j,i]
    }else{
      us.cum[j,i]=t.sims.us[j,i]+us.cum[(j-1),i]
    }
  }
}
mean.us=matrix(NA, 1, 23)
for(i in 1:23){
  mean.us[,i] = mean(us.cum[i,])
}
plot(mean.us[1,], ylim=c(-.25,0))
se.us=matrix(NA,23,1)
for(i in 1:23){
  se.us[i,] = sd(us.cum[i,])
}
upr.us=as.vector(mean.us[1,] + 1.96*se.us[,1])
lines(upr.us,col="lightgrey",type="l")
lwr.us=as.vector(mean.us[1,] - 1.96*se.us[,1])
lines(lwr.us,col="lightgrey",type="l")

all.countries.means[,1]=mean.us[1,]
all.countries.upr[,1]=upr.us
all.countries.lwr[,1]=lwr.us

#
#TUR
#
tur.cum =matrix(NA,23, draws)
t.sims.tur=t(sims.tur)
dim(t.sims.tur)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      tur.cum[j,i] = t.sims.tur[j,i]
    }else{
      tur.cum[j,i]=t.sims.tur[j,i]+tur.cum[(j-1),i]
    }
  }
}
mean.tur=matrix(NA, 1, 23)
for(i in 1:23){
  mean.tur[,i] = mean(tur.cum[i,])
}
plot(mean.tur[1,])
se.tur=matrix(NA,23,1)
for(i in 1:23){
  se.tur[i,] = sd(tur.cum[i,])
}
upr.tur=as.vector(mean.tur[1,] + 1.96*se.tur[,1])
lines(upr.tur,col="lightgrey",type="l")
lwr.tur=as.vector(mean.tur[1,] - 1.96*se.tur[,1])
lines(lwr.tur,col="lightgrey",type="l")

all.countries.means[,18]=mean.tur[1,]
all.countries.upr[,18]=upr.tur
all.countries.lwr[,18]=lwr.tur

#
#IT #14
#
it.cum =matrix(NA,23, draws)
t.sims.it=t(sims.it)
dim(t.sims.it)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      it.cum[j,i] = t.sims.it[j,i]
    }else{
      it.cum[j,i]=t.sims.it[j,i]+it.cum[(j-1),i]
    }
  }
}
mean.it=matrix(NA, 1, 23)
for(i in 1:23){
  mean.it[,i] = mean(it.cum[i,])
}
plot(mean.it[1,])
se.it=matrix(NA,23,1)
for(i in 1:23){
  se.it[i,] = sd(it.cum[i,])
}
upr.it=as.vector(mean.it[1,] + 1.96*se.it[,1])
lines(upr.it,col="lightgrey",type="l")
lwr.it=as.vector(mean.it[1,] - 1.96*se.it[,1])
lines(lwr.it,col="lightgrey",type="l")

all.countries.means[,14]=mean.it[1,]
all.countries.upr[,14]=upr.it
all.countries.lwr[,14]=lwr.it

#
#EL #15
#
el.cum =matrix(NA,23, draws)
t.sims.el=t(sims.el)
dim(t.sims.el)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      el.cum[j,i] = t.sims.el[j,i]
    }else{
      el.cum[j,i]=t.sims.el[j,i]+el.cum[(j-1),i]
    }
  }
}
mean.el=matrix(NA, 1, 23)
for(i in 1:23){
  mean.el[,i] = mean(el.cum[i,])
}
plot(mean.el[1,], ylim=c(-.6,0))
se.el=matrix(NA,23,1)
for(i in 1:23){
  se.el[i,] = sd(el.cum[i,])
}
upr.el=as.vector(mean.el[1,] + 1.96*se.el[,1])
lines(upr.el,col="lightgrey",type="l")
lwr.el=as.vector(mean.el[1,] - 1.96*se.el[,1])
lines(lwr.el,col="lightgrey",type="l")

all.countries.means[,15]=mean.el[1,]
all.countries.upr[,15]=upr.el
all.countries.lwr[,15]=lwr.el

#sims.others1, sims.others2
sims.no = sims.others2[,1:23] #NO #16
dim(sims.no)
sims.dk = sims.others2[,24:46] #DK #17
sims.can = sims.others1[,1:23] #CAN #2
sims.uk = sims.others1[,24:46] #UK #3
sims.nl = sims.others1[,47:69] #NL #4
dim(sims.nl)
sims.be = sims.others1[,70:92] #"Belgium" #[5]   
sims.lu = sims.others1[,93:115] #"Luxembourg" #[6]  
sims.fr = sims.others1[,116:138] #"France" #[7]   
sims.es = sims.others1[,139:161] #"Spain" # [8]  
sims.pt = sims.others1[,162:184] # "Portugal" #[9]   
sims.de = sims.others1[,185:207] # "Germany" #[10]  
sims.pl = sims.others1[,208:230] #"Poland" # [11] 
sims.hu = sims.others1[,231:253] # "Hungary" #[12]  
sims.cz= sims.others1[,254:276] # "Czech Republic" #[13]


#
#NO #16
#
no.cum =matrix(NA,23, draws)
t.sims.no=t(sims.no)
dim(t.sims.no)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      no.cum[j,i] = t.sims.no[j,i]
    }else{
      no.cum[j,i]=t.sims.no[j,i]+no.cum[(j-1),i]
    }
  }
}
mean.no=matrix(NA, 1, 23)
for(i in 1:23){
  mean.no[,i] = mean(no.cum[i,])
}
plot(mean.no[1,])
se.no=matrix(NA,23,1)
for(i in 1:23){
  se.no[i,] = sd(no.cum[i,])
}
upr.no=as.vector(mean.no[1,] + 1.96*se.no[,1])
lines(upr.no,col="lightgrey",type="l")
lwr.no=as.vector(mean.no[1,] - 1.96*se.no[,1])
lines(lwr.no,col="lightgrey",type="l")

all.countries.means[,16]=mean.no[1,]
all.countries.upr[,16]=upr.no
all.countries.lwr[,16]=lwr.no

#
#DK #17
#
dk.cum =matrix(NA,23, draws)
t.sims.dk=t(sims.dk)
dim(t.sims.dk)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      dk.cum[j,i] = t.sims.dk[j,i]
    }else{
      dk.cum[j,i]=t.sims.dk[j,i]+dk.cum[(j-1),i]
    }
  }
}
mean.dk=matrix(NA, 1, 23)
for(i in 1:23){
  mean.dk[,i] = mean(dk.cum[i,])
}
plot(mean.dk[1,])
se.dk=matrix(NA,23,1)
for(i in 1:23){
  se.dk[i,] = sd(dk.cum[i,])
}
upr.dk=as.vector(mean.dk[1,] + 1.96*se.dk[,1])
lines(upr.dk,col="lightgrey",type="l")
lwr.dk=as.vector(mean.dk[1,] - 1.96*se.dk[,1])
lines(lwr.dk,col="lightgrey",type="l")

all.countries.means[,17]=mean.dk[1,]
all.countries.upr[,17]=upr.dk
all.countries.lwr[,17]=lwr.dk

#
#CAN #2
#
can.cum =matrix(NA,23, draws)
t.sims.can=t(sims.can)
dim(t.sims.can)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      can.cum[j,i] = t.sims.can[j,i]
    }else{
      can.cum[j,i]=t.sims.can[j,i]+can.cum[(j-1),i]
    }
  }
}
mean.can=matrix(NA, 1, 23)
for(i in 1:23){
  mean.can[,i] = mean(can.cum[i,])
}
plot(mean.can[1,])
se.can=matrix(NA,23,1)
for(i in 1:23){
  se.can[i,] = sd(can.cum[i,])
}
upr.can=as.vector(mean.can[1,] + 1.96*se.can[,1])
lines(upr.can,col="lightgrey",type="l")
lwr.can=as.vector(mean.can[1,] - 1.96*se.can[,1])
lines(lwr.can,col="lightgrey",type="l")

all.countries.means[,2]=mean.can[1,]
all.countries.upr[,2]=upr.can
all.countries.lwr[,2]=lwr.can

#
#UK #3
#
uk.cum =matrix(NA,23, draws)
t.sims.uk=t(sims.uk)
dim(t.sims.uk)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      uk.cum[j,i] = t.sims.uk[j,i]
    }else{
      uk.cum[j,i]=t.sims.uk[j,i]+uk.cum[(j-1),i]
    }
  }
}
mean.uk=matrix(NA, 1, 23)
for(i in 1:23){
  mean.uk[,i] = mean(uk.cum[i,])
}
plot(mean.uk[1,])
se.uk=matrix(NA,23,1)
for(i in 1:23){
  se.uk[i,] = sd(uk.cum[i,])
}
upr.uk=as.vector(mean.uk[1,] + 1.96*se.uk[,1])
lines(upr.uk,col="lightgrey",type="l")
lwr.uk=as.vector(mean.uk[1,] - 1.96*se.uk[,1])
lines(lwr.uk,col="lightgrey",type="l")

all.countries.means[,3]=mean.uk[1,]
all.countries.upr[,3]=upr.uk
all.countries.lwr[,3]=lwr.uk

#
#NL #4
#
nl.cum =matrix(NA,23, draws)
t.sims.nl=t(sims.nl)
dim(t.sims.nl)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      nl.cum[j,i] = t.sims.nl[j,i]
    }else{
      nl.cum[j,i]=t.sims.nl[j,i]+nl.cum[(j-1),i]
    }
  }
}
mean.nl=matrix(NA, 1, 23)
for(i in 1:23){
  mean.nl[,i] = mean(nl.cum[i,])
}
plot(mean.nl[1,])
se.nl=matrix(NA,23,1)
for(i in 1:23){
  se.nl[i,] = sd(nl.cum[i,])
}
upr.nl=as.vector(mean.nl[1,] + 1.96*se.nl[,1])
lines(upr.nl,col="lightgrey",type="l")
lwr.nl=as.vector(mean.nl[1,] - 1.96*se.nl[,1])
lines(lwr.nl,col="lightgrey",type="l")

all.countries.means[,4]=mean.nl[1,]
all.countries.upr[,4]=upr.nl
all.countries.lwr[,4]=lwr.nl

#
#BE #5
#
be.cum =matrix(NA,23, draws)
t.sims.be=t(sims.be)
dim(t.sims.be)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      be.cum[j,i] = t.sims.be[j,i]
    }else{
      be.cum[j,i]=t.sims.be[j,i]+be.cum[(j-1),i]
    }
  }
}
mean.be=matrix(NA, 1, 23)
for(i in 1:23){
  mean.be[,i] = mean(be.cum[i,])
}
plot(mean.be[1,], ylim=c(-1.5,0))
se.be=matrix(NA,23,1)
for(i in 1:23){
  se.be[i,] = sd(be.cum[i,])
}
upr.be=as.vector(mean.be[1,] + 1.96*se.be[,1])
lines(upr.be,col="lightgrey",type="l")
lwr.be=as.vector(mean.be[1,] - 1.96*se.be[,1])
lines(lwr.be,col="lightgrey",type="l")

all.countries.means[,5]=mean.be[1,]
all.countries.upr[,5]=upr.be
all.countries.lwr[,5]=lwr.be

#
#LU #6
#
lu.cum =matrix(NA,23, draws)
t.sims.lu=t(sims.lu)
dim(t.sims.lu)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      lu.cum[j,i] = t.sims.lu[j,i]
    }else{
      lu.cum[j,i]=t.sims.lu[j,i]+lu.cum[(j-1),i]
    }
  }
}
mean.lu=matrix(NA, 1, 23)
for(i in 1:23){
  mean.lu[,i] = mean(lu.cum[i,])
}
plot(mean.lu[1,])
se.lu=matrix(NA,23,1)
for(i in 1:23){
  se.lu[i,] = sd(lu.cum[i,])
}
upr.lu=as.vector(mean.lu[1,] + 1.96*se.lu[,1])
lines(upr.lu,col="lightgrey",type="l")
lwr.lu=as.vector(mean.lu[1,] - 1.96*se.lu[,1])
lines(lwr.lu,col="lightgrey",type="l")

all.countries.means[,6]=mean.lu[1,]
all.countries.upr[,6]=upr.lu
all.countries.lwr[,6]=lwr.lu

#
#FR #7
#
fr.cum =matrix(NA,23, draws)
t.sims.fr=t(sims.fr)
dim(t.sims.fr)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      fr.cum[j,i] = t.sims.fr[j,i]
    }else{
      fr.cum[j,i]=t.sims.fr[j,i]+fr.cum[(j-1),i]
    }
  }
}
mean.fr=matrix(NA, 1, 23)
for(i in 1:23){
  mean.fr[,i] = mean(fr.cum[i,])
}
plot(mean.fr[1,])
se.fr=matrix(NA,23,1)
for(i in 1:23){
  se.fr[i,] = sd(fr.cum[i,])
}
upr.fr=as.vector(mean.fr[1,] + 1.96*se.fr[,1])
lines(upr.fr,col="lightgrey",type="l")
lwr.fr=as.vector(mean.fr[1,] - 1.96*se.fr[,1])
lines(lwr.fr,col="lightgrey",type="l")

all.countries.means[,7]=mean.fr[1,]
all.countries.upr[,7]=upr.fr
all.countries.lwr[,7]=lwr.fr

#
#ES #8
#
es.cum =matrix(NA,23, draws)
t.sims.es=t(sims.es)
dim(t.sims.es)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      es.cum[j,i] = t.sims.es[j,i]
    }else{
      es.cum[j,i]=t.sims.es[j,i]+es.cum[(j-1),i]
    }
  }
}
mean.es=matrix(NA, 1, 23)
for(i in 1:23){
  mean.es[,i] = mean(es.cum[i,])
}
plot(mean.es[1,])
se.es=matrix(NA,23,1)
for(i in 1:23){
  se.es[i,] = sd(es.cum[i,])
}
upr.es=as.vector(mean.es[1,] + 1.96*se.es[,1])
lines(upr.es,col="lightgrey",type="l")
lwr.es=as.vector(mean.es[1,] - 1.96*se.es[,1])
lines(lwr.es,col="lightgrey",type="l")

all.countries.means[,8]=mean.es[1,]
all.countries.upr[,8]=upr.es
all.countries.lwr[,8]=lwr.es

#
#PT #9
#
pt.cum =matrix(NA,23, draws)
t.sims.pt=t(sims.pt)
dim(t.sims.pt)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      pt.cum[j,i] = t.sims.pt[j,i]
    }else{
      pt.cum[j,i]=t.sims.pt[j,i]+pt.cum[(j-1),i]
    }
  }
}
mean.pt=matrix(NA, 1, 23)
for(i in 1:23){
  mean.pt[,i] = mean(pt.cum[i,])
}
plot(mean.pt[1,])
se.pt=matrix(NA,23,1)
for(i in 1:23){
  se.pt[i,] = sd(pt.cum[i,])
}
upr.pt=as.vector(mean.pt[1,] + 1.96*se.pt[,1])
lines(upr.pt,col="lightgrey",type="l")
lwr.pt=as.vector(mean.pt[1,] - 1.96*se.pt[,1])
lines(lwr.pt,col="lightgrey",type="l")

all.countries.means[,9]=mean.pt[1,]
all.countries.upr[,9]=upr.pt
all.countries.lwr[,9]=lwr.pt

#
#DE #10
#
de.cum =matrix(NA,23, draws)
t.sims.de=t(sims.de)
dim(t.sims.de)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      de.cum[j,i] = t.sims.de[j,i]
    }else{
      de.cum[j,i]=t.sims.de[j,i]+de.cum[(j-1),i]
    }
  }
}
mean.de=matrix(NA, 1, 23)
for(i in 1:23){
  mean.de[,i] = mean(de.cum[i,])
}
plot(mean.de[1,])
se.de=matrix(NA,23,1)
for(i in 1:23){
  se.de[i,] = sd(de.cum[i,])
}
upr.de=as.vector(mean.de[1,] + 1.96*se.de[,1])
lines(upr.de,col="lightgrey",type="l")
lwr.de=as.vector(mean.de[1,] - 1.96*se.de[,1])
lines(lwr.de,col="lightgrey",type="l")

all.countries.means[,10]=mean.de[1,]
all.countries.upr[,10]=upr.de
all.countries.lwr[,10]=lwr.de

#
#PL #11
#
pl.cum =matrix(NA,23, draws)
t.sims.pl=t(sims.pl)
dim(t.sims.pl)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      pl.cum[j,i] = t.sims.pl[j,i]
    }else{
      pl.cum[j,i]=t.sims.pl[j,i]+pl.cum[(j-1),i]
    }
  }
}
mean.pl=matrix(NA, 1, 23)
for(i in 1:23){
  mean.pl[,i] = mean(pl.cum[i,])
}
plot(mean.pl[1,])
se.pl=matrix(NA,23,1)
for(i in 1:23){
  se.pl[i,] = sd(pl.cum[i,])
}
upr.pl=as.vector(mean.pl[1,] + 1.96*se.pl[,1])
lines(upr.pl,col="lightgrey",type="l")
lwr.pl=as.vector(mean.pl[1,] - 1.96*se.pl[,1])
lines(lwr.pl,col="lightgrey",type="l")

all.countries.means[,11]=mean.pl[1,]
all.countries.upr[,11]=upr.pl
all.countries.lwr[,11]=lwr.pl

#
#HU #12
#
hu.cum =matrix(NA,23, draws)
t.sims.hu=t(sims.hu)
dim(t.sims.hu)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      hu.cum[j,i] = t.sims.hu[j,i]
    }else{
      hu.cum[j,i]=t.sims.hu[j,i]+hu.cum[(j-1),i]
    }
  }
}
mean.hu=matrix(NA, 1, 23)
for(i in 1:23){
  mean.hu[,i] = mean(hu.cum[i,])
}
plot(mean.hu[1,])
se.hu=matrix(NA,23,1)
for(i in 1:23){
  se.hu[i,] = sd(hu.cum[i,])
}
upr.hu=as.vector(mean.hu[1,] + 1.96*se.hu[,1])
lines(upr.hu,col="lightgrey",type="l")
lwr.hu=as.vector(mean.hu[1,] - 1.96*se.hu[,1])
lines(lwr.hu,col="lightgrey",type="l")

all.countries.means[,12]=mean.hu[1,]
all.countries.upr[,12]=upr.hu
all.countries.lwr[,12]=lwr.hu

#
#CZ #13
#
cz.cum =matrix(NA,23, draws)
t.sims.cz=t(sims.cz)
dim(t.sims.cz)
for(i in 1:draws){
  for(j in 1:23){
    if(j==1){
      cz.cum[j,i] = t.sims.cz[j,i]
    }else{
      cz.cum[j,i]=t.sims.cz[j,i]+cz.cum[(j-1),i]
    }
  }
}
mean.cz=matrix(NA, 1, 23)
for(i in 1:23){
  mean.cz[,i] = mean(cz.cum[i,])
}
plot(mean.cz[1,])
se.cz=matrix(NA,23,1)
for(i in 1:23){
  se.cz[i,] = sd(cz.cum[i,])
}
upr.cz=as.vector(mean.cz[1,] + 1.96*se.cz[,1])
lines(upr.cz,col="lightgrey",type="l")
lwr.cz=as.vector(mean.cz[1,] - 1.96*se.cz[,1])
lines(lwr.cz,col="lightgrey",type="l")

all.countries.means[,13]=mean.cz[1,]
all.countries.upr[,13]=upr.cz
all.countries.lwr[,13]=lwr.cz

#
#
average=matrix(NA, 23,1)
for(i in 1:23){
  average[i,]=mean(all.countries.means[i,])
}

upr = matrix(NA, 23,1)
for(i in 1:23){
  upr[i,]=mean(all.countries.upr[i,])
}

lwr = matrix(NA, 23,1)
for(i in 1:23){
  lwr[i,]=mean(all.countries.lwr[i,])
}

plot(average[,1], ylim=c(-0.65,0), type="l", lwd=6, xaxt='n', yaxt='n')
lines(upr[,1], col="lightgrey")
lines(lwr[,1], col="lightgrey")

years.list=c(seq(1999,2021,1))
axis(1, at=seq(1,23,1),labels=years.list )

x=as.vector(c(seq(1999,2021,1), seq(2021,1999,-1)))
y=as.vector(c(upr[,1],rev(lwr[,1])))
polygon(x, y, col = "lightgrey")

years=years.list
milsp = upr[,1]
milsp2=lwr[,1]
data=as.data.frame(cbind(years,milsp,milsp2))
data=cbind(data,average)

data0=data
##
##
##
#projection of 15-year period 
data = subset(data, data$years < 2014) #15 years

fig1 <- ggplot(data, aes(x = years)) +
        geom_line(aes(y = average), color = "blue") +  # Main curve
        geom_ribbon(aes(ymin = milsp, ymax = milsp2), alpha = 0.2) +  # Confidence interval
        labs(title = " ",
        x = "Years", y = "Milsp (in %)") +
      theme_minimal()

fig1 <- fig1 + 
  theme(
    axis.title.x = element_text(size = 14),  # Increase the x-axis title size
    axis.title.y = element_text(size = 14)   # Increase the y-axis title size
  )

fig1 <- fig1 + 
  theme(
    panel.border = element_blank(),               
    plot.background = element_rect(fill = "white", color = NA),  
    panel.background = element_rect(fill = "white", color = NA),
    axis.line = element_blank(),                 
    axis.ticks = element_blank()                
  )

print(fig1) 

ggsave("figure1.png", fig1)


#
#average effect of increase in milsp of a single ally
#
average.single=average/18
upr.single=upr/18
lwr.single=lwr/18

plot(average.single[,1], ylim=c(-0.036,0), type="l", lwd=6, , xaxt='n')
lines(upr.single[,1], col="lightgrey")
lines(lwr.single[,1], col="lightgrey")

years.list=c(seq(1999,2021,1))
axis(1, at=seq(1,23,1),labels=years.list )

x=as.vector(c(seq(1999,2021,1), seq(2021,1999,-1)))
y=as.vector(c(upr.single[,1],rev(lwr.single[,1])))
polygon(x, y, col = "lightgrey")

years=years.list
milsp = upr.single[,1]
milsp2=lwr.single[,1]
data=as.data.frame(cbind(years,milsp,milsp2))
data=cbind(data,average.single)

data = subset(data, data$years < 2014)

fig2 <- ggplot(data, aes(x = years)) +
  geom_line(aes(y = average.single), color = "blue") +  # Main curve
  geom_ribbon(aes(ymin = milsp, ymax = milsp2), alpha = 0.2) +  # Confidence interval
  labs(title = " ",
       x = "Years", y = "Milsp (in %)") +
  theme_minimal()

fig2 <- fig2 + 
  theme(
    axis.title.x = element_text(size = 14),  # Increase the x-axis title size
    axis.title.y = element_text(size = 14)   # Increase the y-axis title size
  )

fig2 <- fig2 + 
  theme(
    panel.border = element_blank(),               
    plot.background = element_rect(fill = "white", color = NA),  
    panel.background = element_rect(fill = "white", color = NA),
    axis.line = element_blank(),                 
    axis.ticks = element_blank()                
  )

print(fig2) 

ggsave("figure2.png", fig2)

